#### setup ####

library(modelsummary)
library(interplot)
library(interflex)
library(broom)
library(ggpubr)
library(grid) 
library(tidyverse) 

set.seed(5)

#### Data  ####

data <- read_csv("raw_data.csv")


data <- data %>%
  mutate(
    rid = 1:n(), #  row id
    # Transform party id based on Q3, Q4, Q5
    party7 = case_when(
      Q3 == 1 ~ 0,
      Q3 == 2 ~ 1/6,
      Q5 == 3 ~ 2/6, 
      Q4 == 1 ~ 1,
      Q4 == 2 ~ 5/6,
      Q5 == 1 ~ 4/6,
      Q5 == 2 ~ 3/6
    ),
    party5 = case_when(
      Q3 == 1 ~ 0,
      (Q3 == 2) | (Q5 == 3) ~ 0.25,
      Q4 == 1 ~ 1,
      (Q4 == 2) | (Q5 == 1) ~ 0.75,
      Q5 == 2 ~ 0.5
    ),
    party3 = case_when(
      (Q3 == 2) | (Q5 == 3) | (Q3 == 1) ~ 0,
      (Q4 == 2) | (Q5 == 1) | (Q4 == 1) ~ 1,
      Q5 == 2 ~ 0.5
    ),
    partystrength = case_when(
      Q5 == 2 ~ 0,
      (Q3 == 2) | (Q5 == 3) | (Q4 == 2) | (Q5 == 1) ~ 0.5,
      (Q3 == 1) | (Q4 == 1) ~ 1
    ),
    ind_cond = case_when(
      ind_cond == "rep" ~ 1,
      ind_cond == "dem" ~ 0,
      TRUE ~ NA_real_ 
    ),
    # Combine estimates for rep and dem
    repunderstand = ifelse(format == 'reps_first', Q9, Q12),
    demunderstand = ifelse(format == 'reps_first', Q10, Q11),
    repcompromise = ifelse(format == 'reps_first', Q15, Q18),
    demcompromise = ifelse(format == 'reps_first', Q16, Q17),
    reprightwrong = ifelse(format == 'reps_first', Q21, Q24),
    demrightwrong = ifelse(format == 'reps_first', Q22, Q23),
    repjustfair = ifelse(format == 'reps_first', Q27, Q30),
    demjustfair = ifelse(format == 'reps_first', Q28, Q29),
    # Combine estimates for in and outparty
    inunderstand = ifelse(format == 'reps_first', ifelse(party3 == 1, Q9, Q10), ifelse(party3 == 1, Q12, Q11)),
    outunderstand = ifelse(format == 'reps_first', ifelse(party3 == 1, Q10, Q9), ifelse(party3 == 1, Q11, Q12)),
    incompromise = ifelse(format == 'reps_first', ifelse(party3 == 1, Q15, Q16), ifelse(party3 == 1, Q18, Q17)),
    outcompromise = ifelse(format == 'reps_first', ifelse(party3 == 1, Q16, Q15), ifelse(party3 == 1, Q17, Q18)),
    inrightwrong = ifelse(format == 'reps_first', ifelse(party3 == 1, Q21, Q22), ifelse(party3 == 1, Q24, Q23)),
    outrightwrong = ifelse(format == 'reps_first', ifelse(party3 == 1, Q22, Q21), ifelse(party3 == 1, Q23, Q24)),
    injustfair = ifelse(format == 'reps_first', ifelse(party3 == 1, Q27, Q28), ifelse(party3 == 1, Q30, Q29)),
    outjustfair = ifelse(format == 'reps_first', ifelse(party3 == 1, Q28, Q27), ifelse(party3 == 1, Q29, Q30)),
    treatment= as.numeric(treatment)
  )

data <- data %>%
  mutate(
    compromise_1 = ifelse(str_detect(FL_34_DO, "^Compromise"), 1, 0),
    justfair_1 = ifelse(str_detect(FL_34_DO, "^just/fair"), 1, 0),
    clearrightwrong_1 = ifelse(str_detect(FL_34_DO, "^clearright/wrong"), 1, 0),
    understand_1 = ifelse(str_detect(FL_34_DO, "^Understand"), 1, 0),
    compromise_2 = ifelse(str_detect(FL_34_DO, "\\|Compromise\\|"), 1, 0),
    justfair_2 = ifelse(str_detect(FL_34_DO, "\\|just/fair\\|"), 1, 0),
    clearrightwrong_2 = ifelse(str_detect(FL_34_DO, "\\|clearright/wrong\\|"), 1, 0),
    understand_2 = ifelse(str_detect(FL_34_DO, "\\|Understand\\|"), 1, 0),
    compromise_3 = ifelse(str_detect(FL_34_DO, "Compromise$"), 1, 0),
    justfair_3 = ifelse(str_detect(FL_34_DO, "just/fair$"), 1, 0),
    clearrightwrong_3 = ifelse(str_detect(FL_34_DO, "clearright/wrong$"), 1, 0),
    understand_3 = ifelse(str_detect(FL_34_DO, "Understand$"), 1, 0),
    # reps first is 1
    format = ifelse(format == "reps_first", 1, 0)
  )


data <- data %>%
  mutate(
    understanddiff = ifelse(as.numeric(party7) == 3/6, 
                            (as.numeric(repunderstand) - as.numeric(demunderstand)), 
                            (as.numeric(outunderstand) - as.numeric(inunderstand))),
    compromisediff = ifelse(as.numeric(party7) == 3/6, 
                            (as.numeric(repcompromise) - as.numeric(demcompromise)),
                            (as.numeric(outcompromise) - as.numeric(incompromise))),
    rightwrongdiff = ifelse(as.numeric(party7) == 3/6, 
                            (as.numeric(reprightwrong) - as.numeric(demrightwrong)),
                            (as.numeric(outrightwrong) - as.numeric(inrightwrong))),
    justfairdiff = ifelse(as.numeric(party7) == 3/6, 
                          (as.numeric(repjustfair) - as.numeric(demjustfair)),
                          (as.numeric(outjustfair) - as.numeric(injustfair)))
  )

# outcomes
data <- data %>% 
  mutate(
    inaff = case_when(
      party3 == 1 ~ as.numeric(Q38),
      party3 == 0 ~ as.numeric(Q37),
      TRUE ~ NA_real_
    ),
    outaff = case_when(
      party3 == 1 ~ as.numeric(Q37),
      party3 == 0 ~ as.numeric(Q38),
      TRUE ~ NA_real_
    ),
    netaff = inaff - outaff,
    # Rescale variables so 1 is more frequent discussion/follow government
    # Q7 1-9
    discfreq = case_when(
      Q7 == 9 ~ 1,
      Q7 == 8 ~ 7/8,
      Q7 == 7 ~ 6/8,
      Q7 == 6 ~ 5/8,
      Q7 == 5 ~ 4/8,
      Q7 == 4 ~ 3/8,
      Q7 == 3 ~ 2/8,
      Q7 == 2 ~ 1/8,
      Q7 == 1 ~ 0,
      TRUE ~ NA_real_
    ),
    # Q6 1-4
    followgov = case_when(
      Q6 == 1 ~ 1,
      Q6 == 2 ~ 2/3,
      Q6 == 3 ~ 1/3,
      Q6 == 4 ~ 0,
      TRUE ~ NA_real_
    ), 
    treatchar = ifelse(treatment == 1, "Treatment", "Control")
  )


#### Descriptive #####

comb_understand <- ggplot(data %>% filter(!is.na(party7)), aes(x = factor(party7))) +
  stat_summary(aes(y = repunderstand, color = "repunderstand"), fun = mean, geom = "point", position = position_nudge(x=.2), size =3.5) +
  stat_summary(aes(y = repunderstand, color = "repunderstand"), fun.data = mean_cl_normal, geom = "linerange", position = position_nudge(x = .2), fun.args = list(conf.int = 0.84), size = 1.7) +
  stat_summary(aes(y = repunderstand, color = "repunderstand"), fun.data = mean_cl_normal, geom = "linerange", position = position_nudge(x = .2), fun.args = list(conf.int = 0.95), size = 0.9) +
  stat_summary(aes(y = demunderstand, color = "demunderstand"), fun = mean, geom = "point", size =3.5) +
  stat_summary(aes(y = demunderstand, color = "demunderstand"), fun.data = mean_cl_normal, geom = "linerange", fun.args = list(conf.int = 0.84), size = 1.7) +
  stat_summary(aes(y = demunderstand, color = "demunderstand"), fun.data = mean_cl_normal, geom = "linerange", fun.args = list(conf.int = 0.95), size = 0.9) +
  labs(title = "Asking Me To Understand is an Insult", x = "Respondent Party ID", y = "Mean Estimate (w/ 84 and 95% CIs)") +
  scale_x_discrete(labels = c("Strong Dem", "Weak Dem","Lean Dem","Independent","Lean Rep","Weak Rep", "Strong Rep"), guide = guide_axis(n.dodge = 2)) +
  scale_color_manual(values = c("repunderstand" = "red", "demunderstand" = "blue"))+
  theme_minimal()+
  coord_cartesian(ylim=c(42, 75)) +
  guides(color = FALSE) +
  theme(text = element_text(size=16), plot.title = element_text(face = "italic"))

comb_compromise <- ggplot(data %>% filter(!is.na(party7)), aes(x = factor(party7))) +
  stat_summary(aes(y = repcompromise, color = "repcompromise"), fun = mean, geom = "point", position = position_nudge(x=.2), size =3.5) +
  stat_summary(aes(y = repcompromise, color = "repcompromise"), fun.data = mean_cl_normal, geom = "linerange", position = position_nudge(x = .2), fun.args = list(conf.int = 0.84), size = 1.7) +
  stat_summary(aes(y = repcompromise, color = "repcompromise"), fun.data = mean_cl_normal, geom = "linerange", position = position_nudge(x = .2), fun.args = list(conf.int = 0.95), size = 0.9) +
  stat_summary(aes(y = demcompromise, color = "demcompromise"), fun = mean, geom = "point", size =3.5) +
  stat_summary(aes(y = demcompromise, color = "demcompromise"), fun.data = mean_cl_normal, geom = "linerange", fun.args = list(conf.int = 0.84), size = 1.7) +
  stat_summary(aes(y = demcompromise, color = "demcompromise"), fun.data = mean_cl_normal, geom = "linerange", fun.args = list(conf.int = 0.95), size = 0.9) +
  labs(title = "Compromise Prevents a Better Society", x = "Respondent Party ID", y = "Mean Estimate (w/ 84 and 95% CIs)") +
  scale_x_discrete(labels = c("Strong Dem", "Weak Dem","Lean Dem","Independent","Lean Rep","Weak Rep", "Strong Rep"), guide = guide_axis(n.dodge = 2)) +
  scale_color_manual(values = c("repcompromise" = "red", "demcompromise" = "blue"))+
  theme_minimal()+
  coord_cartesian(ylim=c(42, 75))+
  guides(color = FALSE) +
  theme(text = element_text(size=16), plot.title = element_text(face = "italic"))

comb_justfair <- ggplot(data %>% filter(!is.na(party7)), aes(x = factor(party7))) +
  stat_summary(aes(y = repjustfair, color = "repjustfair"), fun = mean, geom = "point", position = position_nudge(x=.2), size =3.5) +
  stat_summary(aes(y = repjustfair, color = "repjustfair"), fun.data = mean_cl_normal, geom = "linerange", position = position_nudge(x = .2), fun.args = list(conf.int = 0.84), size = 1.7) +
  stat_summary(aes(y = repjustfair, color = "repjustfair"), fun.data = mean_cl_normal, geom = "linerange", position = position_nudge(x = .2), fun.args = list(conf.int = 0.95), size = 0.9) +
  stat_summary(aes(y = demjustfair, color = "demjustfair"), fun = mean, geom = "point", size =3.5) +
  stat_summary(aes(y = demjustfair, color = "demjustfair"), fun.data = mean_cl_normal, geom = "linerange", fun.args = list(conf.int = 0.84), size = 1.7) +
  stat_summary(aes(y = demjustfair, color = "demjustfair"), fun.data = mean_cl_normal, geom = "linerange", fun.args = list(conf.int = 0.95), size = 0.9) +
  labs(title = "My Ideas Are More Just/Fair", x = "Respondent Party ID", y = "Mean Estimate (w/ 84 and 95% CIs)") +
  scale_x_discrete(labels = c("Strong Dem", "Weak Dem","Lean Dem","Independent","Lean Rep","Weak Rep", "Strong Rep"), guide = guide_axis(n.dodge = 2)) +
  scale_color_manual(values = c("repjustfair" = "red", "demjustfair" = "blue"))+
  theme_minimal()+
  coord_cartesian(ylim=c(42, 75))+
  guides(color = FALSE) +
  theme(text = element_text(size=16), plot.title = element_text(face = "italic"))

comb_rightwrong <- ggplot(data %>% filter(!is.na(party7)), aes(x = factor(party7))) +
  stat_summary(aes(y = reprightwrong, color = "reprightwrong"), fun = mean, geom = "point", position = position_nudge(x=.2), size =3.5) +
  stat_summary(aes(y = reprightwrong, color = "reprightwrong"), fun.data = mean_cl_normal, geom = "linerange", position = position_nudge(x = .2), fun.args = list(conf.int = 0.84), size = 1.7) +
  stat_summary(aes(y = reprightwrong, color = "reprightwrong"), fun.data = mean_cl_normal, geom = "linerange", position = position_nudge(x = .2), fun.args = list(conf.int = 0.95), size = 0.9) +
  stat_summary(aes(y = demrightwrong, color = "demrightwrong"), fun = mean, geom = "point", size =3.5) +
  stat_summary(aes(y = demrightwrong, color = "demrightwrong"), fun.data = mean_cl_normal, geom = "linerange", fun.args = list(conf.int = 0.84), size = 1.7) +
  stat_summary(aes(y = demrightwrong, color = "demrightwrong"), fun.data = mean_cl_normal, geom = "linerange", fun.args = list(conf.int = 0.95), size = 0.9) +
  labs(title = "Clear Right and Wrong", x = "Respondent Party ID", y = "Mean Estimate (w/ 84 and 95% CIs)") +
  scale_x_discrete(labels = c("Strong Dem", "Weak Dem","Lean Dem","Independent","Lean Rep","Weak Rep", "Strong Rep"), guide = guide_axis(n.dodge = 2)) +
  scale_color_manual(values = c("reprightwrong" = "red", "demrightwrong" = "blue"))+
  theme_minimal()+
  coord_cartesian(ylim=c(42, 75))+
  guides(color = FALSE) +
  theme(text = element_text(size=16), plot.title = element_text(face = "italic"))

ggarrange(comb_compromise, comb_understand, comb_rightwrong, comb_justfair, ncol = 2, nrow = 2)

# average across attitudes

data %>%
  mutate(
    rep_avg = rowMeans(select(., repcompromise, repunderstand, reprightwrong, repjustfair), na.rm = TRUE),
    dem_avg = rowMeans(select(., demcompromise, demunderstand, demrightwrong, demjustfair), na.rm = TRUE),
    
    party_group = case_when(
      near(party7, 0) | near(party7, 1/6) | near(party7, 2/6) ~ "Democrat",
      near(party7, 1) | near(party7, 5/6) | near(party7, 4/6) ~ "Republican",
      TRUE ~ "Independent"
    )
  ) %>%
  filter(party_group %in% c("Democrat", "Republican")) %>%
  group_by(party_group) %>%
  summarise(
    mean_rep_estimate = mean(rep_avg, na.rm = TRUE),
    sd_rep_estimate = sd(rep_avg, na.rm = TRUE),
    mean_dem_estimate = mean(dem_avg, na.rm = TRUE),
    sd_dem_estimate = sd(dem_avg, na.rm = TRUE)
  )




#### Treatment Analysis ####
data <- data %>%
  mutate(
    dc_dis = (as.numeric(Q35) - 1) / 4,  # standardizing
    dc_agr = (as.numeric(Q36) - 1) / 4,
    dc_diff = dc_agr-dc_dis,
    opdiv = (as.numeric(Q33)-1)/4, 
    hhv = (as.numeric(Q34)-1)/4, 
    inaff_std = inaff/100, 
    outaff_std = outaff/100,
    netaff_std = inaff_std-outaff_std)

data <-  data %>% filter(!is.na(treatchar))

##### ATE Discussion Comfort ######

discdisagree <- ggplot(data = data, aes(x = treatchar, y = as.numeric(dc_dis))) +
  stat_summary(fun.data = mean_cl_normal, geom = "linerange", fun.args = list(conf.int = 0.84), linewidth = 3) +
  stat_summary(fun.data = mean_cl_normal, geom = "linerange", fun.args = list(conf.int = 0.95), linewidth = 1.5) +
  stat_summary(fun = "mean", geom = "point", col = "white", size = 5, shape = 15) +
  stat_summary(aes(label = round(..y.., 3)), fun = mean, geom = "text", size = 4) +
  labs(title = "Discussion Comfort (Disagree)", y = "Mean (w/ 84% and 95% CI)") +
  theme_bw() +
  theme(axis.title.x = element_blank(), axis.text.x = element_text(size = 12), axis.text.y = element_text(size = 10), plot.title = element_text(face = "bold"))

agree_means <- mean(data$dc_agr[data$treatchar == "Treatment"], na.rm = TRUE) -
  mean(data$dc_agr[data$treatchar == "Control"], na.rm = TRUE)

discagree <- ggplot(data = data, aes(x = treatchar, y = as.numeric(dc_agr))) +
  stat_summary(fun.data = mean_cl_normal, geom = "linerange", fun.args = list(conf.int = 0.84), linewidth = 3) +
  stat_summary(fun.data = mean_cl_normal, geom = "linerange", fun.args = list(conf.int = 0.95), linewidth = 1.5) +
  stat_summary(fun = "mean", geom = "point", col = "white", size = 5, shape = 15) +
  stat_summary(aes(label = round(..y.., 3)), fun = mean, geom = "text", size = 4) +
  labs(title = "Discussion Comfort (Agree)", y = "Mean (w/ 84% and 95% CI)") +
  theme_bw() +
  theme(axis.title.x = element_blank(), axis.text.x = element_text(size = 12), axis.text.y = element_text(size = 10), plot.title = element_text(face = "bold"))

discdiff_means <- mean(data$dc_diff[data$treatchar == "Treatment"], na.rm = TRUE) - mean(data$dc_diff[data$treatchar == "Control"], na.rm = TRUE)

discdiff <- ggplot(data = data, aes(x = treatchar, y = dc_diff)) +
  stat_summary(fun.data = mean_cl_normal, geom = "linerange", fun.args = list(conf.int = 0.84), linewidth = 3) +
  stat_summary(fun.data = mean_cl_normal, geom = "linerange", fun.args = list(conf.int = 0.95), linewidth = 1.5) +
  stat_summary(fun = "mean", geom = "point", col = "white", size = 5, shape = 15) +
  stat_summary(aes(label = round(..y.., 3)), fun = mean, geom = "text", size = 4) +
  labs(title = "Likeminded Discussion Preference", y = "Mean (w/ 84% and 95% CI)") +
  theme_bw() +
  theme(axis.title.x = element_blank(), axis.text.x = element_text(size = 12), axis.text.y = element_text(size = 10), plot.title = element_text(face = "bold"))


ATE_discdiff_ttest <- t.test(dc_diff ~ treatchar, data)
ATE_aff_ttest <- t.test(netaff_std ~ treatchar, data)
ATE_hhv_ttest <- t.test(hhv ~ treatchar, data) 
ATE_opdiv_ttest <- t.test(opdiv ~ treatchar, data) 

ATE_df <- data.frame(
  name = c("Discussion Preference", "Affective Polarization", "Helpful Hear Views", "Value Opinion Diversity"),
  point = c(
    ATE_discdiff_ttest$estimate[["mean in group Treatment"]] - ATE_discdiff_ttest$estimate[["mean in group Control"]],
    ATE_aff_ttest$estimate[["mean in group Treatment"]] - ATE_aff_ttest$estimate[["mean in group Control"]],
    ATE_hhv_ttest$estimate[["mean in group Treatment"]] - ATE_hhv_ttest$estimate[["mean in group Control"]], # Using ATE_hhv_ttest (opdiv column)
    ATE_opdiv_ttest$estimate[["mean in group Treatment"]] - ATE_opdiv_ttest$estimate[["mean in group Control"]]  # Using ATE_opdiv_ttest (hhv column)
  ),
  lower_ci = c(
    -ATE_discdiff_ttest$conf.int[2], 
    -ATE_aff_ttest$conf.int[2],
    -ATE_hhv_ttest$conf.int[2],
    -ATE_opdiv_ttest$conf.int[2]
  ),
  upper_ci = c(
    -ATE_discdiff_ttest$conf.int[1], 
    -ATE_aff_ttest$conf.int[1],
    -ATE_hhv_ttest$conf.int[1],
    -ATE_opdiv_ttest$conf.int[1]
  )
)


ATE_disc_plot <- ATE_df %>%
  filter(name == "Discussion Preference") %>%
  ggplot(aes(x = name, y = as.numeric(point))) +
  geom_linerange(aes(ymin = as.numeric(lower_ci), ymax = as.numeric(upper_ci)), linewidth = 1.5) + 
  labs(title = "ATE: Likeminded Discussion Preference", y = "ATE (w/ 95% CI)") +
  theme_bw() +
  theme(axis.title.x = element_blank(), axis.text.x = element_text(size = 12), axis.text.y = element_text(size = 10), plot.title = element_text(face = "bold")) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  stat_summary(fun = "mean", geom = "point", col = "white", size = 5, shape = 15) + 
  geom_point(aes(y=point), colour="white", size=5, shape=15) + 
  geom_text(aes(label = round(as.numeric(point), 3), y=point), size = 4) 

ggarrange(discdiff, ATE_disc_plot, discagree, discdisagree, ncol = 2, nrow = 2)

##### ATE Affect ######



t.test(data$netaff_std[data$treatchar == "Treatment" & !is.na(data$netaff_std)],
       data$netaff_std[data$treatchar == "Control" & !is.na(data$netaff_std)])
t.test(data$inaff_std[data$treatchar == "Treatment" & !is.na(data$inaff_std)],
       data$inaff_std[data$treatchar == "Control" & !is.na(data$inaff_std)])
t.test(data$outaff_std[data$treatchar == "Treatment" & !is.na(data$outaff_std)],
data$outaff_std[data$treatchar == "Control" & !is.na(data$outaff_std)])

inaff_std_plot <- ggplot(data = data, aes(x = treatchar, y = inaff_std)) + 
  stat_summary(fun.data = mean_cl_normal, geom = "linerange", fun.args = list(conf.int = 0.84), linewidth = 3) +
  stat_summary(fun.data = mean_cl_normal, geom = "linerange", fun.args = list(conf.int = 0.95), linewidth = 1.5) +
  stat_summary(fun = "mean", geom = "point", col = "white", size = 5, shape = 15) +
  stat_summary(aes(label = round(..y.., 3)), fun = mean, geom = "text", size = 4) +
  labs(title = "Inparty Feeling Thermometer", y = "Mean (w/ 84% and 95% CI)") +
  theme_bw() +
  theme(axis.title.x = element_blank(), axis.text.x = element_text(size = 12), axis.text.y = element_text(size = 10), plot.title = element_text(face = "bold"))


outaff_std_plot <- ggplot(data = data, aes(x = treatchar, y = outaff_std)) + 
  stat_summary(fun.data = mean_cl_normal, geom = "linerange", fun.args = list(conf.int = 0.84), linewidth = 3) +
  stat_summary(fun.data = mean_cl_normal, geom = "linerange", fun.args = list(conf.int = 0.95), linewidth = 1.5) +
  stat_summary(fun = "mean", geom = "point", col = "white", size = 5, shape = 15) +
  stat_summary(aes(label = round(..y.., 3)), fun = mean, geom = "text", size = 4) +
  labs(title = "Outparty Feeling Thermometer", y = "Mean (w/ 84% and 95% CI)") +
  theme_bw() +
  theme(axis.title.x = element_blank(), axis.text.x = element_text(size = 12), axis.text.y = element_text(size = 10), plot.title = element_text(face = "bold"))


netaff_std_plot <- ggplot(data = data, aes(x = treatchar, y = netaff_std)) + 
  stat_summary(fun.data = mean_cl_normal, geom = "linerange", fun.args = list(conf.int = 0.84), linewidth = 3) +
  stat_summary(fun.data = mean_cl_normal, geom = "linerange", fun.args = list(conf.int = 0.95), linewidth = 1.5) +
  stat_summary(fun = "mean", geom = "point", col = "white", size = 5, shape = 15) +
  stat_summary(aes(label = round(..y.., 3)), fun = mean, geom = "text", size = 4) +
  labs(title = "Inparty - Outparty Affect", y = "Mean (w/ 84% and 95% CI)") +
  theme_bw() +
  theme(axis.title.x = element_blank(), axis.text.x = element_text(size = 12), axis.text.y = element_text(size = 10), plot.title = element_text(face = "bold"))

ATE_aff_plot <- ATE_df %>%
  filter(name == "Affective Polarization") %>%
  ggplot(aes(x = name, y = as.numeric(point))) +
  geom_linerange(aes(ymin = as.numeric(lower_ci), ymax = as.numeric(upper_ci)), linewidth = 1.5) + 
  labs(title = "ATE: Affective Polarization", y = "ATE (w/ 95% CI)") +
  theme_bw() +
  theme(axis.title.x = element_blank(), axis.text.x = element_text(size = 12), axis.text.y = element_text(size = 10), plot.title = element_text(face = "bold")) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_point(aes(y=point), colour="white", size=4, shape=15) + 
  geom_text(aes(label = round(as.numeric(point), 3), y=point), size = 4)

ggarrange(netaff_std_plot, ATE_aff_plot, inaff_std_plot, outaff_std_plot, ncol = 2, nrow = 2)


##### ATE Opinion Diversity ######

t.test(as.numeric(data$opdiv)[data$treatchar == "Treatment" & !is.na(as.numeric(data$opdiv))],
       as.numeric(data$opdiv)[data$treatchar == "Control" & !is.na(as.numeric(data$opdiv))])

opiniondiversity_plot <- ggplot(data = data, aes(x = treatchar, y = opdiv)) + 
  stat_summary(fun.data = mean_cl_normal, geom = "linerange", fun.args = list(conf.int = 0.84), linewidth = 3) +
  stat_summary(fun.data = mean_cl_normal, geom = "linerange", fun.args = list(conf.int = 0.95), linewidth = 1.5) +
  stat_summary(fun = "mean", geom = "point", col = "white", size = 5, shape = 15) +
  stat_summary(aes(label = round(..y.., 3)), fun = mean, geom = "text", size = 4) +
  labs(title = "Diversity of Opinions is Valuable", y = "Mean (w/ 84% and 95% CI)") + 
  theme_bw() +
  theme(axis.title.x = element_blank(), axis.text.x = element_text(size = 12), axis.text.y = element_text(size = 10), plot.title = element_text(face = "bold"))


helphearviews_plot <- ggplot(data = data, aes(x = treatchar, y = hhv)) + 
  stat_summary(fun.data = mean_cl_normal, geom = "linerange", fun.args = list(conf.int = 0.84), linewidth = 3) +
  stat_summary(fun.data = mean_cl_normal, geom = "linerange", fun.args = list(conf.int = 0.95), linewidth = 1.5) +
  stat_summary(fun = "mean", geom = "point", col = "white", size = 5, shape = 15) +
  stat_summary(aes(label = round(..y.., 3)), fun = mean, geom = "text", size = 4) +
  labs(title = "Helpful Hearing Discordant Views", y = "Mean (w/ 84% and 95% CI)") + 
  theme_bw() +
  theme(axis.title.x = element_blank(), axis.text.x = element_text(size = 12), axis.text.y = element_text(size = 10), plot.title = element_text(face = "bold"))

ATE_hhv_plot <- ATE_df %>% 
  filter(name == "Helpful Hear Views") %>%
  ggplot(aes(x = name, y = as.numeric(point))) +
  geom_linerange(aes(ymin = as.numeric(lower_ci), ymax = as.numeric(upper_ci)), linewidth = 1.5) +
  labs(title = "ATE: Helpful Hearing Discordant Views", y = "ATE (w/ 95% CI)") +
  theme_bw() +
  theme(axis.title.x = element_blank(), axis.text.x = element_text(size = 12), axis.text.y = element_text(size = 10), plot.title = element_text(face = "bold")) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_point(aes(y=point), colour="white", size=5, shape=15) +
  geom_text(aes(label = round(as.numeric(point), 3), y=point), size = 4)

ATE_opdiv_plot <- ATE_df %>% 
  filter(name == "Value Opinion Diversity") %>%
  ggplot(aes(x = name, y = as.numeric(point))) +
  geom_linerange(aes(ymin = as.numeric(lower_ci), ymax = as.numeric(upper_ci)), linewidth = 1.5) +
  labs(title = "ATE: Diversity of Opinions is Valuable", y = "ATE (w/ 95% CI)") +
  theme_bw() +
  theme(axis.title.x = element_blank(), axis.text.x = element_text(size = 12), axis.text.y = element_text(size = 10), plot.title = element_text(face = "bold")) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_point(aes(y=point), colour="white", size=5, shape=15) +
  geom_text(aes(label = round(as.numeric(point), 3), y=point), size = 4)

ggarrange(helphearviews_plot, ATE_hhv_plot, opiniondiversity_plot, ATE_opdiv_plot, ncol = 2, nrow = 2)


#### Heterogeneity (followgov) ####

###### Affect DV (followgov) ######

summary(netaff_std_int_fg <- lm(netaff_std ~ treatment + followgov + treatment * followgov, data = data))

netaff_stdreg_fg_plot <- modelplot(netaff_std_int_fg, coef_rename = c("treatment" = "Treatment", "followgov" = "Follow Gov", "treatment:followgov" = "Treatment * Follow Gov")) +
  labs(x = 'Coefficients', y = 'Term names', title = 'Affective Polarization (Follow Gov)') +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold"), plot.title.position = 'plot') +
  geom_vline(xintercept = 0, linetype = 2)

netaff_stdINT_fg_plot <- interplot(m = netaff_std_int_fg, var1 = "treatment", var2 = "followgov", hist = TRUE) +
  labs(x = 'Follow Gov/Public Affairs', y = 'Marginal Effect of Treatment', title = 'Affective Polarization') +
  theme_bw() +
  theme(plot.title = element_text(face = "bold"), axis.title = element_text(size = 12)) +
  geom_hline(yintercept = 0, linetype = "dashed")

###### Discussion DV (followgov) ######

summary(disdif_int_fg <- lm(as.numeric(dc_diff) ~ treatment + followgov + treatment * followgov, data = data))

discdiffreg_fg_plot <- modelplot(disdif_int_fg, coef_rename = c("treatment" = "Treatment", "followgov" = "Follow Gov", "treatment:followgov" = "Treatment * Follow Gov")) +
  labs(x = 'Coefficients', y = 'Term names', title = 'Likeminded Discussion Preference (Follow Gov)') +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold"), plot.title.position = 'plot') +
  geom_vline(xintercept = 0, linetype = 2)

discdiffINT_fg_plot <- interplot(m = disdif_int_fg, var1 = "treatment", var2 = "followgov", hist = TRUE) +
  labs(x = 'Follow Gov/Public Affairs', y = 'Marginal Effect of Treatment', title = 'Likeminded Discussion Preference') +
  theme_bw() +
  theme(plot.title = element_text(face = "bold"), axis.title = element_text(size = 12)) +
  geom_hline(yintercept = 0, linetype = "dashed")

###### Opinion Diversity DV ######

summary(opdiv_int_fg <- lm(opdiv ~ treatment + followgov + treatment * followgov, data = data)) 

opiniondiversityreg_fg_plot <- modelplot(opdiv_int_fg, coef_rename = c("treatment" = "Treatment", "followgov" = "Follow Gov", "treatment:followgov" = "Treatment * Follow Gov")) +
  labs(x = 'Coefficients', y = 'Term names', title = 'Value Opinion Diversity (Follow Gov)') +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold"), plot.title.position = 'plot') +
  geom_vline(xintercept = 0, linetype = 2)

opiniondiversityINT_fg_plot <- interplot(m = opdiv_int_fg, var1 = "treatment", var2 = "followgov", hist = TRUE) +
  labs(x = 'Follow Gov/Public Affairs', y = 'Marginal Effect of Treatment', title = 'Value Opinion Diversity') +
  theme_bw() +
  theme(plot.title = element_text(face = "bold"), axis.title = element_text(size = 12)) +
  geom_hline(yintercept = 0, linetype = "dashed")

summary(hhv_int_fg <- lm(hhv ~ treatment + followgov + treatment * followgov, data = data)) 

helphearviewsreg_fg_plot <- modelplot(hhv_int_fg, coef_rename = c("treatment" = "Treatment", "followgov" = "Follow Gov", "treatment:followgov" = "Treatment * Follow Gov")) +
  labs(x = 'Coefficients', y = 'Term names', title = 'Helpful Hearing Discordant Views (Follow Gov)') +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold"), plot.title.position = 'plot') +
  geom_vline(xintercept = 0, linetype = 2)

helphearviewsINT_fg_plot <- interplot(m = hhv_int_fg, var1 = "treatment", var2 = "followgov", hist = TRUE) +
  labs(x = 'Follow Gov/Public Affairs', y = 'Marginal Effect of Treatment', title = 'Helpful Hearing Discordant Views') +
  theme_bw() +
  theme(plot.title = element_text(face = "bold"), axis.title = element_text(size = 12)) +
  geom_hline(yintercept = 0, linetype = "dashed")


#### Heterogeneity (discussion frequency) #####


###### Affect DV (discussion frequency) #####

summary(netaff_std_int_df <- lm(netaff_std ~ treatment + discfreq + treatment * discfreq, data = data))

Dnetaff_stdreg_df_plot <- modelplot(netaff_std_int_df, coef_rename = c("treatment" = "Treatment", "discfreq" = "Discussion Freq", "treatment:discfreq" = "Treatment * Discussion Freq")) +
  labs(x = 'Coefficients', y = 'Term names', title = 'Affective Polarization (Discussion Freq)') +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold"), plot.title.position = 'plot') +
  geom_vline(xintercept = 0, linetype = 2)

Dnetaff_stdINT_df_plot <- interplot(m = netaff_std_int_df, var1 = "treatment", var2 = "discfreq", hist = TRUE) +
  labs(x = 'Discussion Frequency', y = 'Marginal Effect of Treatment', title = 'Affective Polarization') +
  theme_bw() +
  theme(plot.title = element_text(face = "bold"), axis.title = element_text(size = 12)) +
  geom_hline(yintercept = 0, linetype = "dashed")

###### Discussion DV (discussion frequency) ######

summary(disdif_int_df <- lm(as.numeric(dc_diff) ~ treatment + discfreq + treatment * discfreq, data = data)) 

Ddiscdiffreg_df_plot <- modelplot(disdif_int_df, coef_rename = c("treatment" = "Treatment", "discfreq" = "Discussion Freq", "treatment:discfreq" = "Treatment * Discussion Freq")) +
  labs(x = 'Coefficients', y = 'Term names', title = 'Likeminded Discussion Preference (Discussion Freq)') +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold"), plot.title.position = 'plot') +
  geom_vline(xintercept = 0, linetype = 2)

DdiscdiffINT_df_plot <- interplot(m = disdif_int_df, var1 = "treatment", var2 = "discfreq", hist = TRUE) +
  labs(x = 'Discussion Frequency', y = 'Marginal Effect of Treatment', title = 'Likeminded Discussion Preference') +
  theme_bw() +
  theme(plot.title = element_text(face = "bold"), axis.title = element_text(size = 12)) +
  geom_hline(yintercept = 0, linetype = "dashed")

###### Opinion Diversity DV (discussion frequency) ######
summary(opdiv_int_df <- lm(opdiv ~ treatment + discfreq + treatment * discfreq, data = data))

Dopiniondiversityreg_df_plot <- modelplot(opdiv_int_df, coef_rename = c("treatment" = "Treatment", "discfreq" = "Discussion Freq", "treatment:discfreq" = "Treatment * Discussion Freq")) +
  labs(x = 'Coefficients', y = 'Term names', title = 'Value Opinion Diversity (Discussion Freq)') +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold"), plot.title.position = 'plot') +
  geom_vline(xintercept = 0, linetype = 2)

DopiniondiversityINT_df_plot <- interplot(m = opdiv_int_df, var1 = "treatment", var2 = "discfreq", hist = TRUE) +
  labs(x = 'Discussion Frequency', y = 'Marginal Effect of Treatment', title = 'Value Opinion Diversity') +
  theme_bw() +
  theme(plot.title = element_text(face = "bold"), axis.title = element_text(size = 12)) +
  geom_hline(yintercept = 0, linetype = "dashed")

summary(hhv_int_df <- lm(hhv ~ treatment + discfreq + treatment * discfreq, data = data)) 

Dhelphearviewsreg_df_plot <- modelplot(hhv_int_df, coef_rename = c("treatment" = "Treatment", "discfreq" = "Discussion Freq", "treatment:discfreq" = "Treatment * Discussion Freq")) +
  labs(x = 'Coefficients', y = 'Term names', title = 'Helpful Hearing Discordant Views (Discussion Freq)') +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold"), plot.title.position = 'plot') +
  geom_vline(xintercept = 0, linetype = 2)

DhelphearviewsINT_df_plot <- interplot(m = hhv_int_df, var1 = "treatment", var2 = "discfreq", hist = TRUE) +
  labs(x = 'Discussion Frequency', y = 'Marginal Effect of Treatment', title = 'Helpful Hearing Discordant Views') +
  theme_bw() +
  theme(plot.title = element_text(face = "bold"), axis.title = element_text(size = 12)) +
  geom_hline(yintercept = 0, linetype = "dashed")

#### Heterogeneity Summaries #####

affmodels_list <- list("Follow Gov Heterogeneity" = netaff_std_int_fg, "Discussion Freq. Heterogeneity" = netaff_std_int_df)
modelsummary(affmodels_list, coef_map = list("treatment" = "Treatment",
                                             "followgov" = "Follow Gov.",
                                             "treatment:followgov" = "Treatment X Follow Gov.",
                                             "discfreq" = "Discussion Freq",
                                             "treatment:discfreq" = "Treatment X Discussion Freq"),
             gof_omit = "AIC|BIC|Log|F|RMSE", estimate = "{estimate}{stars}", output = "latex")


discmodels_list <- list("Follow Gov Heterogeneity" = disdif_int_fg, "Discussion Freq. Heterogeneity" = disdif_int_df)
modelsummary(discmodels_list, coef_map = list("treatment" = "Treatment",
                                              "followgov" = "Follow Gov.",
                                              "treatment:followgov" = "Treatment X Follow Gov.",
                                              "discfreq" = "Discussion Freq",
                                              "treatment:discfreq" = "Treatment X Discussion Freq"),
             gof_omit = "AIC|BIC|Log|F|RMSE", estimate = "{estimate}{stars}", output = "latex")


# opdiv (Q33) and hhv (Q34) model summaries
opdivmodels_list <- list(
  "Value Opinion Diversity (Follow Gov)" = opdiv_int_fg, 
  "Value Opinion Diversity (Disc Freq)" = opdiv_int_df  
)
modelsummary(opdivmodels_list, coef_map = list("treatment" = "Treatment",
                                               "followgov" = "Follow Gov.", "treatment:followgov" = "Treatment X Follow Gov.",
                                               "discfreq" = "Discussion Freq", "treatment:discfreq" = "Treatment X Discussion Freq"),
             gof_omit = "AIC|BIC|Log|F|RMSE", estimate = "{estimate}{stars}", output = "latex")


hhvmodels_list <- list(
  "Helpful Hearing Views (Follow Gov)" = hhv_int_fg, 
  "Helpful Hearing Views (Disc Freq)" = hhv_int_df  
)
modelsummary(hhvmodels_list, coef_map = list("treatment" = "Treatment",
                                             "followgov" = "Follow Gov.", "treatment:followgov" = "Treatment X Follow Gov.",
                                             "discfreq" = "Discussion Freq", "treatment:discfreq" = "Treatment X Discussion Freq"),
             gof_omit = "AIC|BIC|Log|F|RMSE", estimate = "{estimate}{stars}", output = "latex")


models_int <- list(
  "Affective Polarization (Follow Gov.)" = netaff_std_int_fg,
  "Affective Polarization (Talk Freq.)" = netaff_std_int_df, 
  "Discussion Preference (Follow Gov.)" = disdif_int_fg,
  "Discussion Preference (Talk Freq.)" = disdif_int_df,  
  "Value Opinion Diversity (Follow Gov.)" = opdiv_int_fg, 
  "Value Opinion Diversity (Talk Freq.)" = opdiv_int_df, 
  "Helpful Hearing Views (Follow Gov.)" = hhv_int_fg,   
  "Helpful Hearing Views (Talk Freq.)" = hhv_int_df)

model_order <- names(models_int)
tidy_models <- lapply(models_int, function(model) {
  tidy(model, conf.int = TRUE)
})

for (i in seq_along(tidy_models)) {
  tidy_models[[i]]$model <- names(models_int)[i]
}

coef_df <- bind_rows(tidy_models) %>% filter(term != "(Intercept)")

coef_df <- coef_df %>%
  mutate(term = recode(term,
                       "treatment" = "Treatment",
                       "discfreq" = "Political Talk Frequency",
                       "followgov" = "Typically Follow Government",
                       "treatment:discfreq" = "Treatment x Talk Freq.",
                       "treatment:followgov" = "Treatment x Follow Gov."))

coef_order <- c("Political Talk Frequency", "Typically Follow Government",
                "Treatment", "Treatment x Talk Freq.", "Treatment x Follow Gov.")

coef_df <- coef_df %>%
  mutate(model = factor(model, levels = model_order),
         term = factor(term, levels = coef_order))

ggplot(coef_df, aes(x = term, y = estimate,
                    ymin = conf.low, ymax = conf.high,
                    color = model)) +
  geom_point(position = position_dodge2(width = 1, padding = 0.2), size = 5) + 
  geom_linerange(position = position_dodge2(width = 1, padding = 0.2), linewidth = 2) + 
  scale_color_viridis_d(option = "C", end = 0.9, begin = 0.05)+
  coord_flip() +
  labs(title = "Treatment Heterogeneity Summary", x = "", y = "", color = "Key: DV (Treatment Moderator)") +
  theme_bw() +
  theme(
    text = element_text(size = 16),
    axis.text = element_text(size = 16),
    axis.title = element_text(size = 16),
    plot.title = element_text(size = 18, face = "bold"),
    legend.title = element_text(size = 18, face = "bold"),
    legend.text = element_text(size = 16)
  ) +
  geom_hline(yintercept = 0, linetype = "dashed")


# individual het plots 

# Affective Polarization
bothnetaff_std <- ggarrange(netaff_stdINT_fg_plot, netaff_stdreg_fg_plot, Dnetaff_stdINT_df_plot, Dnetaff_stdreg_df_plot, ncol = 2, nrow = 2)

# Opinion Diversity (Value Opinion Diversity - opdiv/Q33)
bothopdiv <- ggarrange(opiniondiversityINT_fg_plot, opiniondiversityreg_fg_plot, DopiniondiversityINT_df_plot, Dopiniondiversityreg_df_plot, ncol = 2, nrow = 2)

# Helpful Hearing Views (hhv/Q34)
bothhhv <- ggarrange(helphearviewsINT_fg_plot, helphearviewsreg_fg_plot, DhelphearviewsINT_df_plot, Dhelphearviewsreg_df_plot, ncol = 2, nrow = 2)

# Discussion Difference
bothdiscdiff <- ggarrange(discdiffINT_fg_plot, discdiffreg_fg_plot, DdiscdiffINT_df_plot, Ddiscdiffreg_df_plot, ncol = 2, nrow = 2)


#### Learning Analysis  #### 

data <- data %>%
  mutate(
    est_in_1 = case_when(
      compromise_1 == 1 ~ as.numeric(incompromise),
      understand_1 == 1 ~ as.numeric(inunderstand),
      justfair_1 == 1 ~ as.numeric(injustfair),
      clearrightwrong_1 == 1 ~ as.numeric(inrightwrong),
      TRUE ~ NA_real_ 
    ),
    est_out_1 = case_when(
      compromise_1 == 1 ~ as.numeric(outcompromise),
      understand_1 == 1 ~ as.numeric(outunderstand),
      justfair_1 == 1 ~ as.numeric(outjustfair),
      clearrightwrong_1 == 1 ~ as.numeric(outrightwrong),
      TRUE ~ NA_real_
    ),
    est_in_2 = case_when(
      compromise_2 == 1 ~ as.numeric(incompromise),
      understand_2 == 1 ~ as.numeric(inunderstand),
      justfair_2 == 1 ~ as.numeric(injustfair),
      clearrightwrong_2 == 1 ~ as.numeric(inrightwrong),
      TRUE ~ NA_real_
    ),
    est_out_2 = case_when(
      compromise_2 == 1 ~ as.numeric(outcompromise),
      understand_2 == 1 ~ as.numeric(outunderstand),
      justfair_2 == 1 ~ as.numeric(outjustfair),
      clearrightwrong_2 == 1 ~ as.numeric(outrightwrong),
      TRUE ~ NA_real_
    ),
    est_in_3 = case_when(
      compromise_3 == 1 ~ as.numeric(incompromise),
      understand_3 == 1 ~ as.numeric(inunderstand),
      justfair_3 == 1 ~ as.numeric(injustfair),
      clearrightwrong_3 == 1 ~ as.numeric(inrightwrong),
      TRUE ~ NA_real_
    ),
    est_out_3 = case_when(
      compromise_3 == 1 ~ as.numeric(outcompromise),
      understand_3 == 1 ~ as.numeric(outunderstand),
      justfair_3 == 1 ~ as.numeric(outjustfair),
      clearrightwrong_3 == 1 ~ as.numeric(outrightwrong),
      TRUE ~ NA_real_
    ),
    first_est_cat = case_when(
      compromise_1 == 1 ~ "compromise",
      understand_1 == 1 ~ "understand",
      justfair_1 == 1 ~ "justfair",
      clearrightwrong_1 == 1 ~ "rightwrong",
      TRUE ~ NA_character_ 
    ),
    second_est_cat = case_when(
      compromise_2 == 1 ~ "compromise",
      understand_2 == 1 ~ "understand",
      justfair_2 == 1 ~ "justfair",
      clearrightwrong_2 == 1 ~ "rightwrong",
      TRUE ~ NA_character_ 
    ),
    third_est_cat = case_when(
      compromise_3 == 1 ~ "compromise",
      understand_3 == 1 ~ "understand",
      justfair_3 == 1 ~ "justfair",
      clearrightwrong_3 == 1 ~ "rightwrong",
      TRUE ~ NA_character_ 
    ),
    omitted = case_when(
      compromise_1 == 0 & compromise_2 == 0 & compromise_3 == 0 & !is.na(first_est_cat) ~ "compromise",
      understand_1 == 0 & understand_2 == 0 & understand_3 == 0 & !is.na(first_est_cat) ~ "understand",
      justfair_1 == 0 & justfair_2 == 0 & justfair_3 == 0 & !is.na(first_est_cat) ~ "justfair",
      clearrightwrong_1 == 0 & clearrightwrong_2 == 0 & clearrightwrong_3 == 0 & !is.na(first_est_cat) ~ "rightwrong",
      TRUE ~ NA_character_ 
    )
  )



justtreat <- filter(data, treatment==1)
justcontrol <- filter(data, treatment==0)



in_treat <- data.frame(
  x = c("First set", "Second set", "Third set"),
  y = c(mean_cl_normal(justtreat$est_in_1)$y, mean_cl_normal(justtreat$est_in_2)$y, mean_cl_normal(justtreat$est_in_3)$y),  
  ymin = c(mean_cl_normal(justtreat$est_in_1)$ymin, mean_cl_normal(justtreat$est_in_2)$ymin, mean_cl_normal(justtreat$est_in_3)$ymin),  
  ymax = c(mean_cl_normal(justtreat$est_in_1)$ymax, mean_cl_normal(justtreat$est_in_2)$ymax, mean_cl_normal(justtreat$est_in_3)$ymax))

in_control <- data.frame(
  x = c("First set", "Second set", "Third set"),
  y = c(mean_cl_normal(justcontrol$est_in_1)$y, mean_cl_normal(justcontrol$est_in_2)$y, mean_cl_normal(justcontrol$est_in_3)$y),  
  ymin = c(mean_cl_normal(justcontrol$est_in_1)$ymin, mean_cl_normal(justcontrol$est_in_2)$ymin, mean_cl_normal(justcontrol$est_in_3)$ymin),  
  ymax = c(mean_cl_normal(justcontrol$est_in_1)$ymax, mean_cl_normal(justcontrol$est_in_2)$ymax, mean_cl_normal(justcontrol$est_in_3)$ymax))

out_treat <- data.frame(
  x = c("First set", "Second set", "Third set"),
  y = c(mean_cl_normal(justtreat$est_out_1)$y, mean_cl_normal(justtreat$est_out_2)$y, mean_cl_normal(justtreat$est_out_3)$y), 
  ymin = c(mean_cl_normal(justtreat$est_out_1)$ymin, mean_cl_normal(justtreat$est_out_2)$ymin, mean_cl_normal(justtreat$est_out_3)$ymin),  
  ymax = c(mean_cl_normal(justtreat$est_out_1)$ymax, mean_cl_normal(justtreat$est_out_2)$ymax, mean_cl_normal(justtreat$est_out_3)$ymax))

out_control <- data.frame(
  x = c("First set", "Second set", "Third set"),
  y = c(mean_cl_normal(justcontrol$est_out_1)$y, mean_cl_normal(justcontrol$est_out_2)$y, mean_cl_normal(justcontrol$est_out_3)$y),  
  ymin = c(mean_cl_normal(justcontrol$est_out_1)$ymin, mean_cl_normal(justcontrol$est_out_2)$ymin, mean_cl_normal(justcontrol$est_out_3)$ymin),  
  ymax = c(mean_cl_normal(justcontrol$est_out_1)$ymax, mean_cl_normal(justcontrol$est_out_2)$ymax, mean_cl_normal(justcontrol$est_out_3)$ymax))


learn_1 <- list(t.test(justtreat$est_in_1, justcontrol$est_in_1)) 
learn_1_vector <- unlist(learn_1) # Convert list to vector
learn_2 <- list(t.test(justtreat$est_in_2, justcontrol$est_in_2))
learn_2_vector <- unlist(learn_2) # Convert list to vector
learn_3 <- list(t.test(justtreat$est_in_3, justcontrol$est_in_3))
learn_3_vector <- unlist(learn_3) # Convert list to vector
learn_4 <- list(t.test(justtreat$est_out_1, justcontrol$est_out_1))
learn_4_vector <- unlist(learn_4) # Convert list to vector
learn_5 <- list(t.test(justtreat$est_out_2, justcontrol$est_out_2))
learn_5_vector <- unlist(learn_5) # Convert list to vector
learn_6 <- list(t.test(justtreat$est_out_3, justcontrol$est_out_3))
learn_6_vector <- unlist(learn_6) # Convert list to vector

learn_df <- data.frame(
  set = factor(c("First Estimate", "Second Estimate", "Third Estimate", 
                 "First Estimate", "Second Estimate", "Third Estimate"), 
               levels = c("First Estimate", "Second Estimate", "Third Estimate")), # Ensure order for lines
  target = c("In", "In", "In", "Out", "Out", "Out"),
  point = c(
    as.numeric(learn_1[[1]][["estimate"]][["mean of x"]] - learn_1[[1]][["estimate"]][["mean of y"]]),
    as.numeric(learn_2[[1]][["estimate"]][["mean of x"]] - learn_2[[1]][["estimate"]][["mean of y"]]),
    as.numeric(learn_3[[1]][["estimate"]][["mean of x"]] - learn_3[[1]][["estimate"]][["mean of y"]]),
    as.numeric(learn_4[[1]][["estimate"]][["mean of x"]] - learn_4[[1]][["estimate"]][["mean of y"]]),
    as.numeric(learn_5[[1]][["estimate"]][["mean of x"]] - learn_5[[1]][["estimate"]][["mean of y"]]),
    as.numeric(learn_6[[1]][["estimate"]][["mean of x"]] - learn_6[[1]][["estimate"]][["mean of y"]])
  ),
  lower_ci = as.numeric(c(learn_1_vector[4], learn_2_vector[4], learn_3_vector[4], learn_4_vector[4], learn_5_vector[4], learn_6_vector[4])),
  upper_ci = as.numeric(c(learn_1_vector[5], learn_2_vector[5], learn_3_vector[5], learn_4_vector[5], learn_5_vector[5], learn_6_vector[5]))
)

ggplot() +
  geom_point(data = learn_df, aes(x = set, y = point, color = target), size = 4, position = position_dodge2(width = 0.1)) +
  geom_linerange(data = learn_df, aes(x = set, y = point, ymin = lower_ci, ymax = upper_ci, color = target), position = position_dodge2(width = 0.1)) +
  geom_line(data = filter(learn_df, target == "In"), aes(x = set, y = point, group = target, color = target), position = position_dodge2(width = 0.1)) +
  geom_line(data = filter(learn_df, target == "Out"), aes(x = set, y = point, group = target, color = target), position = position_dodge2(width = 0.1)) +
  scale_color_manual(values = c("In" = "#009E73", "Out" = "darkorange2")) +
  labs(
    x = "Set",
    y = "Difference in Estimates Between Control and Treatment",
    color = "Estimated Party",
    title = "Increasing Divergence Between Treatment and Control Estimates"
  ) +
  theme_bw() +
  theme(text = element_text(size = 12)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "darkgrey")



#### Observational Extension ####

# standardizes 0-1
data <- data %>%
  filter(!is.na(est_in_1) & !is.na(est_out_1)) %>% 
  mutate(
    est_in_1 = est_in_1 / 100, 
    est_out_1 = est_out_1 / 100,
    first_est_cat = as.factor(first_est_cat)
  )

coef_names <- c(
  "est_out_1" = "First Out-group Estimate (0-1)", 
  "est_in_1" = "First In-group Estimate (0-1)",  
  "treatment" = "Treatment",           
  "first_est_catrightwrong" = "Right/Wrong (vs. Compromise)", 
  "first_est_catjustfair" = "Just/Fair (vs. Compromise)",
  "first_est_catunderstand" = "Understand (vs. Compromise)",
  "est_out_1:est_in_1" = "Out-group Est. × In-group Est."
)

###### Discussion Comfort DVs ######
summary(dis_dis_obs <- lm(dc_dis ~ est_out_1 * est_in_1 + treatment + first_est_cat, data = data))
summary(dis_agr_obs <- lm(dc_agr ~ est_out_1 * est_in_1 + treatment + first_est_cat, data = data))
summary(dis_dif_obs <- lm(dc_diff ~ est_out_1 * est_in_1 + treatment + first_est_cat, data = data))

# table
models_dc <- list(
  "DV: Comfort w/ Disagreer" = dis_dis_obs, 
  "DV: Comfort w/ Agreer" = dis_agr_obs,
  "DV: Likeminded Discussion Pref." = dis_dif_obs
)
modelsummary(models_dc, output = "kableExtra", stars = TRUE, coef_map = coef_names, title = "Observational Models: Discussion Comfort Outcomes")

# interflex plots 
data$first_est_cat <- as.factor(data$first_est_cat) 

dis_dis_obs_het_out <- interflex(Y = "dc_dis", D = "est_out_1", X = "est_in_1", Z = c("treatment", "first_est_cat"),
                                 data = as.data.frame(data), estimator = "kernel", nboots = 5000, neval=101, 
                                 parallel = TRUE, na.rm = TRUE, Xlabel = "In-Party Perception (0-1)", theme.bw = TRUE,
                                 main = "DV: Discussion Comfort (Disagree)", cex.axis = 0.8, cex.lab = 0.9,
                                 cex.main = 0.8, ylab = "Marginal Effect of Out-Party Perception")

dis_agr_obs_het_out <- interflex(Y = "dc_agr", D = "est_out_1", X = "est_in_1", Z = c("treatment", "first_est_cat"),
                                 data = as.data.frame(data), estimator = "kernel", nboots = 5000, neval=101,
                                 parallel = TRUE, na.rm = TRUE, Xlabel = "In-Party Perception (0-1)", theme.bw = TRUE,
                                 main = "DV: Discussion Comfort (Agree)", cex.axis = 0.8, cex.lab = 0.9,
                                 cex.main = 0.8, ylab = "Marginal Effect of Out-Party Perception")

dis_dif_obs_het_out <- interflex(Y = "dc_diff", D = "est_out_1", X = "est_in_1", Z = c("treatment", "first_est_cat"),
                                 data = as.data.frame(data), estimator = "kernel", nboots = 5000, neval = 101, 
                                 parallel = TRUE, na.rm = TRUE, Xlabel = "In-Party Perception (0-1)", theme.bw = TRUE,
                                 main = "DV: Likeminded Discussion Preference", cex.axis = 0.8, cex.lab = 0.9,
                                 cex.main = 0.8, ylab = "Marginal Effect of Out-Party Perception")



##### Affect DVs####
summary(aff_out_obs <- lm(outaff_std ~ est_out_1 * est_in_1 + treatment + first_est_cat, data = data))
summary(aff_in_obs  <- lm(inaff_std ~ est_out_1 * est_in_1 + treatment + first_est_cat, data = data))
summary(aff_net_obs <- lm(netaff_std ~ est_out_1 * est_in_1 + treatment + first_est_cat, data = data))

# table
models_aff <- list(
  "DV: Out-Party Affect (0-1)" = aff_out_obs,
  "DV: In-Party Affect (0-1)" = aff_in_obs,
  "DV: Affective Polarization (scaled)" = aff_net_obs
)
modelsummary(models_aff, output = "kableExtra", stars = TRUE, coef_map = coef_names, title = "Observational Models: Affective Outcomes")

# Interflex Plots
aff_out_obs_het_out <- interflex(Y = "outaff_std", D = "est_out_1", X = "est_in_1", Z = c("treatment", "first_est_cat"),
                                 data = as.data.frame(data), estimator = "kernel", nboots = 5000, neval=101,
                                 parallel = TRUE, na.rm = TRUE, Xlabel = "In-Party Perception (0-1)", theme.bw = TRUE,
                                 main = "DV: Out-Party Affect (0-1)", cex.axis = 0.8, cex.lab = 0.9,
                                 cex.main = 0.8, ylab = "Marginal Effect of Out-Party Perception")

aff_in_obs_het_out <- interflex(Y = "inaff_std", D = "est_out_1", X = "est_in_1", Z = c("treatment", "first_est_cat"),
                                data = as.data.frame(data), estimator = "kernel", nboots = 5000, neval=101,
                                parallel = TRUE, na.rm = TRUE, Xlabel = "In-Party Perception (0-1)", theme.bw = TRUE,
                                main = "DV: In-Party Affect (0-1)", cex.axis = 0.8, cex.lab = 0.9,
                                cex.main = 0.8, ylab = "Marginal Effect of Out-Party Perception")

aff_net_obs_het_out <- interflex(Y = "netaff_std", D = "est_out_1", X = "est_in_1", Z = c("treatment", "first_est_cat"),
                                 data = as.data.frame(data), estimator = "kernel", nboots = 5000, neval=101,
                                 parallel = TRUE, na.rm = TRUE, Xlabel = "In-Party Perception (0-1)", theme.bw = TRUE,
                                 main = "DV: Affective Polarization (scaled)", cex.axis = 0.8, cex.lab = 0.9,
                                 cex.main = 0.8, ylab = "Marginal Effect of Out-Party Perception")


##### Opinion Diversity DV #####
summary(hhv_obs   <- lm(hhv ~ est_out_1 * est_in_1 + treatment + first_est_cat, data = data))
summary(opdiv_obs <- lm(opdiv ~ est_out_1 * est_in_1 + treatment + first_est_cat, data = data))

# Table
models_hhdiv <- list(
  "DV: Helpful Hearing Discordant Views (0-1)" = hhv_obs,
  "DV: Value Opinion Diversity (0-1)" = opdiv_obs
)
modelsummary(models_hhdiv, output = "kableExtra", stars = TRUE, coef_map = coef_names, title = "Observational Models: Opinion Diversity Outcomes")

# Interflex Plots
hhv_obs_het_out <- interflex(Y = "hhv", D = "est_out_1", X = "est_in_1", Z = c("treatment", "first_est_cat"),
                             data = as.data.frame(data), estimator = "kernel", nboots = 5000, neval=101,
                             parallel = TRUE, na.rm = TRUE, Xlabel = "In-Party Perception (0-1)", theme.bw = TRUE,
                             main = "DV: Helpful Hearing Discordant Views (0-1)", cex.axis = 0.8, cex.lab = 0.9,
                             cex.main = 0.8, ylab = "Marginal Effect of Out-Party Perception")

opdiv_obs_het_out <- interflex(Y = "opdiv", D = "est_out_1", X = "est_in_1", Z = c("treatment", "first_est_cat"),
                               data = as.data.frame(data), estimator = "kernel", nboots = 5000, neval=101,
                               parallel = TRUE, na.rm = TRUE, Xlabel = "In-Party Perception (0-1)", theme.bw = TRUE,
                               main = "DV: Value Opinion Diversity (0-1)", cex.axis = 0.8, cex.lab = 0.9,
                               cex.main = 0.8, ylab = "Marginal Effect of Out-Party Perception")


##### combined plots ####

# Affect Outcomes 
cb_palette_aff <- c("Affective Polarization (scaled)" = "#6A3D9A", "In-Party Affect (0-1)" = "#009E73", "Out-Party Affect (0-1)" = "#E69F00")

coef_out_aff <- tidy(aff_out_obs) %>% mutate(model = "Out-Party Affect (0-1)") 
coef_in_aff  <- tidy(aff_in_obs)  %>% mutate(model = "In-Party Affect (0-1)")  
coef_net_aff <- tidy(aff_net_obs) %>% mutate(model = "Affective Polarization (scaled)") 

coef_data_aff <- bind_rows(coef_net_aff, coef_in_aff, coef_out_aff) %>%
  filter(!(term == "(Intercept)" | grepl("^first_est_cat", term))) 

coef_data_aff$model <- factor(coef_data_aff$model, levels = c("Out-Party Affect (0-1)", "In-Party Affect (0-1)", "Affective Polarization (scaled)"))

coef_plot_aff <- ggplot(coef_data_aff, aes(x = term, y = estimate, color = model)) +
  geom_pointrange(aes(ymin = estimate - 1.96 * std.error, ymax = estimate + 1.96 * std.error),
  position = position_dodge(width = 0.5), size = 0.4) +
  coord_flip() +
  scale_color_manual(values = cb_palette_aff, guide = "none") +
  scale_x_discrete(labels = c( 
    "treatment" = "Treatment Condition",
    "est_in_1" = "In-party Estimate (0-1)",
    "est_out_1" = "Out-party Estimate (0-1)",
    "est_out_1:est_in_1" = "In-party X Out-party Est."
  )) +
  theme_bw() +
  theme( axis.title.x = element_blank(), 
         axis.title.y = element_blank(),
         axis.text.x = element_text(size = 14), 
         axis.text.y = element_text(size = 14),
         plot.title = element_text(hjust = 0, size = 16, face="bold")) +
  geom_hline(yintercept = 0, linetype = "dotted", size = 0.5) +
  labs(title = "Observational: Affect Outcomes")

# Arrange and print Affect plots
grid.newpage()
pushViewport(viewport(layout = grid.layout(2, 2)))
print(coef_plot_aff,       vp = viewport(layout.pos.row = 1, layout.pos.col = 1))
print(aff_net_obs_het_out, vp = viewport(layout.pos.row = 1, layout.pos.col = 2))
print(aff_in_obs_het_out,  vp = viewport(layout.pos.row = 2, layout.pos.col = 1))
print(aff_out_obs_het_out, vp = viewport(layout.pos.row = 2, layout.pos.col = 2))


# Discussion Comfort Outcomes 
cb_palette_dis <- c("Likeminded Discussion Pref." = "#6A3D9A", "Comfort w/ Agreer" = "#009E73", "Comfort w/ Disagreer" = "#E69F00")

coef_dis_dis <- tidy(dis_dis_obs) %>% mutate(model = "Comfort w/ Disagreer")
coef_dis_agr <- tidy(dis_agr_obs) %>% mutate(model = "Comfort w/ Agreer")
coef_dis_dif <- tidy(dis_dif_obs) %>% mutate(model = "Likeminded Discussion Pref.")

coef_data_dis <- bind_rows(coef_dis_dif, coef_dis_agr, coef_dis_dis) %>% 
  filter(!(term == "(Intercept)" | grepl("^first_est_cat", term)))

coef_data_dis$model <- factor(coef_data_dis$model, levels = c("Comfort w/ Disagreer", "Comfort w/ Agreer", "Likeminded Discussion Pref."))

coef_plot_dis <- ggplot(coef_data_dis, aes(x = term, y = estimate, color = model)) +
  geom_pointrange(aes(ymin = estimate - 1.96 * std.error, ymax = estimate + 1.96 * std.error),
                  position = position_dodge(width = 0.5), size = 0.4) +
  coord_flip() +
  scale_color_manual(values = cb_palette_dis, guide = "none") +
  scale_x_discrete(labels = c(
    "treatment" = "Treatment Condition",
    "est_in_1" = "In-party Estimate (0-1)",
    "est_out_1" = "Out-party Estimate (0-1)",
    "est_out_1:est_in_1" = "In-party X Out-party Est."
  )) +
  theme_bw() +
  theme( axis.title.x = element_blank(), 
         axis.title.y = element_blank(),
         axis.text.x = element_text(size = 14), 
         axis.text.y = element_text(size = 14),
         plot.title = element_text(hjust = 0, size = 16, face="bold")) +
  geom_hline(yintercept = 0, linetype = "dotted", size = 0.5) +
  labs(title = "Observational: Discussion Comfort")

# Arrange and print Discussion Comfort plots
grid.newpage()
pushViewport(viewport(layout = grid.layout(2, 2)))
print(coef_plot_dis,       vp = viewport(layout.pos.row = 1, layout.pos.col = 1))
print(dis_dif_obs_het_out, vp = viewport(layout.pos.row = 1, layout.pos.col = 2))
print(dis_agr_obs_het_out, vp = viewport(layout.pos.row = 2, layout.pos.col = 1))
print(dis_dis_obs_het_out, vp = viewport(layout.pos.row = 2, layout.pos.col = 2))


# Opinion Diversity Outcomes 
cb_palette_hhdv <- c("Helpful Hearing Discordant Views (0-1)" = "#009E73", "Value Opinion Diversity (0-1)" = "#E69F00")

coef_hhv   <- tidy(hhv_obs)   %>% mutate(model = "Helpful Hearing Discordant Views (0-1)")
coef_opdiv <- tidy(opdiv_obs) %>% mutate(model = "Value Opinion Diversity (0-1)")

coef_data_hhdv <- bind_rows(coef_hhv, coef_opdiv) %>%
  filter(!(term == "(Intercept)" | grepl("^first_est_cat", term)))

coef_data_hhdv$model <- factor(coef_data_hhdv$model, levels = c("Helpful Hearing Discordant Views (0-1)", "Value Opinion Diversity (0-1)"))

coef_plot_hhdv <- ggplot(coef_data_hhdv, aes(x = term, y = estimate, color = model)) +
  geom_pointrange(aes(ymin = estimate - 1.96 * std.error, ymax = estimate + 1.96 * std.error),
                  position = position_dodge(width = 0.5), size = 0.4) +
  coord_flip() +
  scale_color_manual(values = cb_palette_hhdv, guide = "none") +
  scale_x_discrete(labels = c(
    "treatment" = "Treatment Condition",
    "est_in_1" = "In-party Estimate (0-1)",
    "est_out_1" = "Out-party Estimate (0-1)",
    "est_out_1:est_in_1" = "In-party X Out-party Est."
  )) +
  theme_bw() +
  theme( axis.title.x = element_blank(), 
         axis.title.y = element_blank(),
         axis.text.x = element_text(size = 14), 
         axis.text.y = element_text(size = 14),
         plot.title = element_text(hjust = 0, size = 16, face="bold")) +
  geom_hline(yintercept = 0, linetype = "dotted", size = 0.5) +
  labs(title = "Observational: Utility in Opposing Ideas")

# Arrange and print Opinion Diversity plots
grid.newpage()
pushViewport(viewport(layout = grid.layout(2, 2)))
print(coef_plot_hhdv,    vp = viewport(layout.pos.row = 1, layout.pos.col = 1))
print(hhv_obs_het_out,   vp = viewport(layout.pos.row = 2, layout.pos.col = 2))
print(opdiv_obs_het_out, vp = viewport(layout.pos.row = 2, layout.pos.col = 1))